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Abstract 

d ' Validation and verification of plasma kinetics codes requires the development of 

Q_i' quantitative methods and techniques for code comparisons. We describe two pa- 

c/J , rameters that can be used for characterization of differences between such codes. 

It is shown that these parameters, which are determined from the most general re- 
£^' suits of kinetic codes, can provide important information on the differences between 

(— I . the basic rate coefficients employed. Application of this method is illustrated by 

OhI comparisons of some results from the 3rd NLTE Code Comparison Workshop for 

carbon, germanium, and gold plasmas. 
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1 Introduction 



Spectroscopy of plasmas has a wide range of applications (|l|; |2) , including the 
study of astrophysical objects, diagnostics of laser produced plasmas, and anal- 

^ \ ysis of radiation emission from plumes of rockets. In spite of the importance 

Fri of this subject and the substantial efforts made by numerous groups, there 

persist significant discrepancies in results of the plasma kinetic codes used to 

^ ' analyze plasma emission spectra under non-local-thermodynamic-equilibrium 

(NLTE) conditions. For examples of the level of disagreement see, e.g., the 
reports from the NLTE Code Comparison Workshops (|3|;y;|5|). 

Given the plasma particle temperature and density, the first task of kinetic 
codes is the computation of the charge state and excitation state distributions 
dJiy). These are typically determined from the set of rate equations: 
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all populating processes 

/ , n^Nf^mMcm^Cm', (1) 

all depopulating processes 

where A''^^ is the density of ions having charge C ( < C ^ ^ with Z being the 
nuclear charge) excited to state m (ordered according to their ascending en- 
ergy, m = corresponds to the ion's ground state), Ue is the electron density, 
and M^a^b denotes the rate coefficient (r-c) for the transition of an ion from 
state a to state b due to some atomic process. In Eq. (1), k represents the 
number of electrons taking part in any given process, thus, k = for sponta- 
neous decay and autoionization, /c = 1 for the electron impact processes (e.g., 
excitation and photorecombination), and k = 2 for the three-body recombina- 
tion. In large-size plasmas, where linear dimension is larger than the photon 
mean free path, photon induced processes must also be included in (1). In a 
steady state plasma, dN(^rn/dt = 0, and Eq. (1) reduces to: 

all populating processes all depopulating processes 

This is a finite set of non-linear coupled equations for N(;m 's. If one is interested 
only in the density of the charge states, regardless of the ionic excitations, Eq. 
(2) can be further simplified, and the result is a set of recursive equations: 

Nr r(2) p(3) • y'^) 



(2) 

Here A*"^ = J2m^(;m is the density of the charge state, and J^^^+i, i^A^i^^ 

and R(-^i^^ are the total rate coefficients for ionization, two-body (radia- 
tive-|-dielectronic) recombination, and three-body recombination, respectively. 
For quasineutral plasmas, solutions of Eq.(3) have to satisfy two complemen- 
tary conditions, namely, 

z z 

ni = J2 N^; ne = Y.CNc = ^^- (4) 

C=o c=o 



Here rii is the total ion density, and Z is the average charge of the ions. In 
the following we assume that N^ is the fractional abundance of charge state 
(, which is equivalent to the assumption that rii = 1. 

The disagreement between the results of the kinetic codes developed by various 
researchers may originate from several factors. One of the most important is 
the approximate character of the r-c's. In the literature one can find several 



recommended formulas for the relevant r-c's. Some may be more accurate than 
others, but none was shown to have high accuracy. For higher accuracy, the r- 
c's can be directly calculated for each transition between atomic states, using 
advanced quantum-mechanical methods. The necessity to generate a large 
number of r-c's for kinetic calculations, however, impedes the application of 
these techniques, thereby forcing a compromise between computational speed 
and accuracy. The solutions of (3) obviously depend on the methods chosen 
for the determination of the r-c's in the right-hand side of (3). 

Another source of disagreement is the criterion used for the continuum or 
ionization potential lowering in plasmas. This is particularly important in 
high-density plasmas where the plasma potential moves the upper ionic bound 
states into the continuum, leaving the ion only with a finite number of discrete 
states. Moreover, due to the fluctuations of the local microfield around each 
ion, even the "tightly bound states" may change into instantaneous quasi- 
molecular states whose treatment is, as yet, not clear. 

Comparison of the results of NLTE kinetic codes was a subject of several 
workshops (|3|;y;|5|). The participants were asked to submit large sets of various 
physical quantities to be compared, and in numerous cases very significant 
differences were found. This situation is well exemplified in Fig. 1, where the 
various calculations of the relative ion populations for a germanium plasma 
are presented for a specific case from the NLTE-3 Workshop (J5[). One can 
clearly notice a significant spread both for the mean ion charges and for the 
distribution widths calculated with different kinetic codes. Although a variety 
of physical parameters investigated at the Workshop was mostly sufficient 
to draw conclusions about sources of discrepancies, introduction of simple 
and clearly defined new parameters that would pinpoint some fundamental 
underlying differences would greatly facilitate such comparisons. 



In Ref. (jTl) one of the authors developed a method that allows determination of 
variations in an ionization state distribution due to small changes in the rates 
of basic atomic processes. The purpose of the present paper is the opposite, 
that is, to introduce new parameters that quantify differences in the input 
atomic rates in kinetic codes using the results of calculations. Such approach 
is similar to solution of the inverse problem in physics. 

Below we address a simple but critical question, namely: "How different are 
two population kinetic codes?" It is obvious that the most comprehensive an- 
swer can be obtained provided all input and output parameters as well as the 
complete description of the approximations used are available. Unfortunately, 
this is not always the case. Here we introduce parameters that are straightfor- 
ward to calculate and can provide an answer in a clear manner. An important 
feature of the proposed method is that only the most general kinetic charac- 
teristics, such as the mean ion charge and central momenta, are required to 
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Fig. 1. Relative ion populations for a steady-state germanium plasma at electron 
temperature Tg = 250 eV and electron density Ng = 10^^ cm~^ (|5l;la)- 

determine these new parameters. 

The paper is organized as foUows. In the next section we define these pa- 
rameters and explain their meanings. In Section 3 we apply the method to 
characterize the differences between code results for three representative cases 
from the NLTE-3 Workshop (|5|). Finally, in Section 4, a short summary and 
conclusions are presented. 



2 Quantitative characterization of kinetic codes 



We introduce the notation Q^ for the r.h.s. of Eq. (3): 



Qc 



'C--C+1 



r(2) , „ p{3) 



(5) 



and denote by AQ^ = Qref — Q the overall difference in this quantity for a 
given code relative to the reference. Using this definition, Eq.(3) is rewritten 



as: 



iVc+i 



N, 



Qc- (6) 



c 



It was shown in Ref. p[), that the relative difference in the final results, 
AN(^/N(^, caused by the relative differences used in the input data, AQ^/Q^, 
is given by 



AA^^ 



A^. 



C 



a^ — a, (7) 



where 

«c = z. -7s— ; « = z. «c— • (8) 

('=0 ^C c=o * 

In particular, if the relative differences are all equal, then 

^Qc/Qc — V = const, (9) 

and Eq.(7) reduces to the simple form (J7|): 

p(C-^). (10) 



AAT^ 



In Ref. (|7|) a full discussion is presented about the meaning of Eq. (10), as well 
as how and when this difference influences the results of kinetic codes. Even if 
AQ^/Q^ is not constant, one can define their average value, p = (AQ^/Q^)^, 

over the relevant ion charge states. 

The above formulas were derived assuming small differences between code 
results. In a general case, however, the deviations are not small, and one has 
to symmetrize the relevant parameters with respect to both compared codes, 
i.e.: 

iVc = ^(^i,C + ^2,c), Z=^(Z, + Z,), (11) 



and so on, which is assumed in what follows. 



Multiplying both sides of Eq.(lO) by {( — Z)N(^ and summing over all charge 
states, one arrives at: 



AZ 

P = = =2 

Z^-Z 



AZ 

0-2 



(12) 



where o"2 is the variance, or second central moment, which is related to the 
ionization distribution width. Remarkably, this equation links the average dif- 
ference in the ratios of effective atomic rates to the most general plasma pa- 
rameters, namely, mean ion charges and variances. 

However, if two compared codes have the same Z but different variances, the 
parameter p defined by Eq. (12) is zero. This simply means that the approx- 
imation (9) is insensitive to differences in ionization distribution widths. To 
take this dependence into account, one can add the next term in the expansion 
with respect to {( — Z): 



AQ^ 

Qc 



p + 2k{C-Z). 



(13) 



Substituting (13) into Eqs. (8) and (7), one obtains: 






piC -Z) + k (C-Zy + iC-Z)- a2 



(14) 



To derive equations relating two parameters, p and k, we multiply Eq. (14) 
by {( — Z)N(^ and ((^ — Z^Ni^ and then sum both sides over ion charges. The 
two ensuing equations are sufficient to obtain the following expressions: 



p 



AZ (0-4 + 0-3 - o-f ) - Ao-2 (0-2 



0-3) 



0-4 (72 

Aa2 • o'2 — AZ ■ o"3 



^i 



al 



0-40-2 



O-q 



0-3 

^2 



(15) 



( — Z) N(- is the ith central moment. 



where a^ = Ef=o' 

Equation (15) is the main result of the present paper. We propose to use 
the parameters p and k for characterization of differences between plasma 
kinetic codes. Using Eq. (15), which depends only on the most general kinetic 
parameters, namely, mean ion charges and central momenta, one can directly 
evaluate the average differences between the effective rates implemented in 
various kinetic codes. 



3 Comparison of kinetic codes 



The above described method is apphed here to the computational results from 
the 3rd NLTE Workshop (J5|) that can be accessed in the NIST SAHA database 
(J6|). This database contains various parameters, including mean ion charges, 
central momenta, and ion populations, which may be used to determine the 
quantities p and k defined in Eqs. (12) and (15). The SAHA database also 
provides other valuable parameters, such as effective rates and partition func- 
tions, so that a user can obtain a deep insight into differences between kinetic 
codes. In accordance with the policy accepted by the Workshop participants, 
the results will be presented without direct attribution, although a list of par- 
ticipating codes will be given for each case. Note also that not all codes provide 
a complete set of central momenta up to 174. 

The first step of the comparison procedure consists in selection of a reference 
against which the other codes are to be compared. In the following comparisons 
the reference code is chosen arbitrarily, as generally there are no a priori 
physical grounds to prefer a particular code. Obviously, the average values of 
p and k would change when selecting another code as a reference. However, 
the standard deviations a of the corresponding distributions of p and k that 
refiect the average spread within a group of codes should not change and 
therefore are reported below as well. In what follows, the parameters p and 
k determined from Eqs. (12) and (15) are referred to as "calculated", while 
those determined by fitting the AQ/Q ratios (Eqs. 9 and 13) are referred to 
as "fit" . As the fitting procedure has to include only the physically significant 
cases, ion states with A^^ < 10"*^ were excluded from the comparison. Also, 
in order to emphasize the contribution from the most populated states, the 
fitting was performed with the weights g^ = -y/iViiN^, where Ni and N2 are 
the ion populations of the reference and compared codes. 

Among numerous cases available in the SAHA database we selected three 
cases for germanium, carbon, and gold. These elements cover a wide range 
of ion charges and their ions represent atomic systems with different level of 
complexity. 



3.1 Ge 



For germanium, we selected a relatively simple case of Tg = 600 eV and 
Ng = 10^^ cm~'^, where almost all codes have a mean ion charge Z ^ 22 
corresponding to a closed-shell Ne-like Ge. This case will be discussed in more 
detail than the C and Au cases. 

Table 1 presents calculated and fit (superscript "/" ) values of parameters p 



Table 1 

Calculated and fit parameters p and k for the Ge case of Tg = 600 eV and Ne = 
10^^ cm~^. Superscript / denotes fit values. Subscript denotes p's determined from 
the one-parameter formulas (9) and (12). Code 7 was excluded in determination of 
standard deviation a. 



Code No. 


Po 


Po 


P 


k 


P^ 


kf 


1 


-0.566 


-0.708 


-0.509 


-0.249 


-0.510 


-0.308 


2 


-0.400 


-0.247 


-0.379 


-0.044 


-0.309 


0.124 


3 


-0.507 


-0.417 


-0.571 


0.130 


-0.594 


0.288 


4 


0.223 


0.222 


0.159 


0.150 


0.140 


0.256 


5 


-0.203 


-0.140 


-0.210 


0.033 


-0.206 


0.113 


6 


-0.764 


-0.942 


-0.760 


-0.130 


-0.760 


-0.260 


7 


-4.266 


-1.298 


-3.800 


-1.389 






8 


-0.592 


-0.122 


-0.334 


-0.175 


-0.018 


0.413 


a 


0.301 


0.362 


0.272 


0.142 


0.298 


0.255 



and k. The values of p determined from the single-parameter formula (12) or 
fit using Eq. (9) have subscript "0". One can immediately notice a generally 
good agreement between the calculated and fit values of p for all but one code. 
Table 1 clearly demonstrates that code 7 is an outlier, which is also emphasized 
by its very different value of the mean ion charge Z k, 27. Moreover, the two- 
parameter fit was not performed for code 7, as for only one ion stage both 
code 7 and the reference code have populations larger than 10"^. 

Another interesting feature is a very small difference between the "simple" p of 
Eq. (9) and calculated and fit values of p determined from the two-parameter 
formulas. While agreement between differently calculated parameters p is gen- 
erally very good, the calculated and fit values of k show worse level of corre- 
spondence. This is not surprising since /c is a high-order parameter which may 
be more sensitive to small variations in data. 

As already mentioned, a standard deviation a would unambiguously represent 
the spread of parameters p and k within a particular group of codes. The last 
row in Table 1 show cr's calculated for each column (the outlier code 7 was 
not included in the determination of a). Remarkably, a's for calculated pq and 
two-parameter calculated and fit p's agree within only 6 %, and even o"(po) 
deviates by less than 25 %. The value of a ^ 0.3 means that in this group of 
codes the average deviation of effective ionization and recombination rates is 
about 30 %. 



Finally, consider the dependence of a on electron temperature Tg. The SAHA 
database contains data for T^ = 150 eV, 250 eV, 450 eV, and 600 eV at Ue 
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Fig. 2. Standard deviation a for calculated parameters po (circles), p (squares), 
and k (triangles) as a function of electron temperature for the germanium case at 
electron density Ne 



IQi'^ cm-3. 



= 10 cm" . The calculated a's for pq, p, and k are presented in Fig. 2 for 
the four temperatures. At low Tg, the mostly populated ions are those with 
open shells, and since these cases are most difficult to calculate, the difference 
between codes is the largest. With the increase of electron temperature, the 
ionization stage approaches the closed-shell Ne-like ion, and therefore agree- 
ment improves dramatically. Note also that a^ is smaller than ap, and therefore 
in many cases a simple one-parameter formula for p would be sufficient for es- 
timates of the difference between codes. 



3.2 C 



Carbon cases in the SAHA database were calculated at a single electron den- 
sity of Ne = 10^^ cm~^. At Tg = 20 eV, which is the selected case here, the 
mean ion charge varies between 1.8 to 3.3 for different codes. There are no 
obvious outliers in this case, and therefore two-parameter fits were successful 
for all codes (Table 2). Similar to the Ge case, the correlation between cal- 
culated and fit values is generally good, although one-parameter p's for codes 
8 and 9 show larger discrepancies. The standard deviations for the different 
parameters p in Table 2 agree even better than for the Ge case, however, they 
are about a factor of 2 larger and reach 60 %. This would seem to be un- 



Table 2 

Calculated and fit parameters p and k for the C case of Tg = 20 eV and Ng = 10^^ 
cin~^. Superscript / denotes fit values. Subscript denotes p's determined from the 
one-parameter formulas (9) and (12). 



Code No. 


Po 


pI 


P 


k 


p' 


kf 


1 


-0.795 


-1.003 


-1.182 


0.411 


-1.440 


0.693 


2 


1.043 


0.537 


0.755 


0.515 


0.635 


0.605 


3 


0.629 


0.434 


0.443 


0.295 


0.435 


0.361 


4 


-1.015 


-0.991 


-1.245 


0.154 


-1.046 


0.135 


5 


-0.770 


-0.899 


-1.064 


0.291 


-1.001 


0.305 


6 


-0.345 


-0.337 


-0.393 


0.042 


-0.340 


0.009 


7 


-0.209 


-0.395 


-0.467 


0.331 


-0.471 


0.374 


8 


0.020 


-0.158 


-0.221 


0.322 


-0.218 


0.403 


9 


0.257 


-0.374 


-0.132 


0.547 


-0.362 


0.851 


a 


0.651 


0.535 


0.660 


0.151 


0.638 


0.250 



expected as carbon simulations could include only seven ionization stages at 
most and are thus supposed to be simpler. The reason for the larger o"'s is that 
the effects of ionization potential lowering are much more important here due 
to the significantly higher density, and different treatments of the continuum 
lowering noticeably contribute to the increased spread of the p and, to a lesser 
extent, k values. 



3.3 Au 



The gold cases available in the SAHA database show very significant differ- 
ences, and for the selected case of Ng = 10^^ cm~^ and Tg = 750 eV, the 
mean ion charge varies between 31 (code 7) and 44 (code 4). Moreover, ioniza- 
tion distributions for codes 1 and 2 show double peak structure unlike other, 
smoother bell-like distributions. It is therefore not surprising that the absolute 
values of p, which are the main indicators of code disagreements, are much 
larger than for the Ge and C cases discussed above. As a consequence, the 
standard deviations for the differently determined p and k values do not show 
the same level of agreement as previously. This situation simply reflects very 
significant differences between code results submitted to the 3rd NLTE Code 
Comparison Workshop. 
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Table 3 

Calculated and fit parameters p and k for the Au case of Tg = 750 eV and Ng = 10^^ 
cm"^. Superscript / denotes fit values. Subscript denotes jo's determined from the 
one-parameter formulas (9) and (12). 



Code No. 


Po 


A 


P 


k 


P^ 


kf 


1 


2.293 


1.376 


1.118 


-0.454 


1.297 


0.042 


2 


-0.534 


-0.410 


0.418 


-0.085 


-0.434 


-0.025 


3 


-1.151 


-1.061 


-0.393 


0.050 


-1.067 


0.050 


4 


-1.338 


-1.178 


-0.497 


0.061 


-1.189 


0.052 


5 


1.638 


1.387 


0.384 


-0.101 


1.419 


-0.039 


6 


-1.086 


-1.042 


-0.422 


0.070 


-1.052 


0.069 


7 


3.393 


-0.970 


0.678 


-0.129 


-0.979 


-0.056 


8 


-1.162 


-0.302 


-0.302 


-0.028 


-0.218 


0.403 


0" 


1.763 


1.001 


0.568 


0.160 


0.996 


0.136 



4 Conclusions 



As the complexity of plasma kinetic codes rapidly increases, their verification 
and validation is becoming mandatory for establishing credibility of compu- 
tational results. To this end, a development of new techniques for code com- 
parisons is an urgent and important task. In the present paper we introduced 
two new parameters for the characterization of discrepancies between plasma 
kinetic codes. These parameters describe differences between effective ioniza- 
tion and recombination rates used in the codes. Importantly, the only physical 
quantities required for their calculation are the mean ion charges and central 
momenta that are the most widely reported characteristics of plasma kinetic 
calculations. Since the final formulas include only the simplest algebra, this 
method provides very fast estimates of code differences in the input atomic 
rates. The new parametrization was applied to the data from the 3rd NLTE 
Code Comparison Workshop and the presented results clearly prove simplicity 
and reliability of the method used. We plan to implement this method to the 
analysis of the data from future NLTE workshops. 
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